Ground state and thermal properties of a lattice gas on a cylindrical surface 
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Adsorbed gases within, or outside of, carbon nanotubes may be analyzed with an approximate 
model of adsorption on lattice sites situated on a cylindrical surface. Using this model, the ground 
state energies of alternative lattice structures are calculated, assuming Lennard- Jones pair interac- 
tions between the particles. The resulting energy and equilibrium structure are nonanalytic functions 
of radius (R) because of commensuration effects associated with the cylindrical geometry. Specifi- 
cally, as R varies, structural transitions occur between configurations differing in the "ring number", 
defined as the number of atoms located at a common value of the longitudinal coordinate (z). The 
thermodynamic behavior of this system is evaluated at finite temperatures, using a Hamiltonian 
with nearest-neighbor interactions. The resulting specific heat bears a qualitative resemblance to 
that of the one-dimensional Ising model. 
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INTRODUCTION 

Atoms or molecules may be confined either within, or 
on the outside surface of, cylindrical materials, such as 
carbon nanotubes. The existence of such systems raises a 
set of interesting questions concerning the thermal, struc- 
tural and dynamical properties of a cylindrical monolayer 
film. In this paper, we evaluate the thermal and struc- 
tural properties of such a film, with the help of two sim- 
plified models. We believe that some of our results are 
realistic, although others may be artifacts of the model 
|1, 1^. For example. Figure 1, taken from previous work 
1 2 1 , depicts the density of H2 molecules within a nanotube 
of radius R=7 A, at temperature T=10 K. Note that the 
radial spread of the so-called "cylindrical shell" phase is 
some 10% of the mean radial distance (r) « 3.8 A. In 
such a situation, the model of confinement to a precise 
value of R would seem appropriate. 

The first task we undertake is to ascertain the ground 
state energy and structure of an ensemble of atoms, 
assumed to be classical, that interact with all other 
atoms with a Lennard- Jones (L-J) pair potential: U{r) = 
4e[((T/r)-^^ — (cr/r)^]. Here a is the nominal diameter of 
the atoms and e is the well-depth of their mutual inter- 
action. We assume, without proof, that the equilibrium 
structure is a periodic crystal. To determine its proper- 
ties, we minimize the ground state energy per particle in 
this cylindrical surface lattice. In so doing, we consider 
possible periodic structures and determine that structure 
which has the lowest energy, at any specified value of R. 
The second problem is to evaluate the thermal properties 
of such a "lattice gas" which has varying fractional site 
occupancy 9. In this case, we simplify the problem by 
including just nearest neighbor interactions. 

The behavior of the system can be expressed in terms 
of a reduced radius R* = R/a. Both the ground state 
and finite temperature (T) problems have much-studied 
one-dimensional (ID) and 2D limits, corresponding to 
R* = and 00, respectively. Interestingly, the behavior 



does not interpolate smoothly between these limits as R* 
is varied. This happens because of a commensuration 
effect arising from the periodicity associated with the az- 
imuthal angle (j). This phenomenon is analogous to that 
found in monolayer films in the regime where the film's 
lattice constant is similar to that of the substrate [Q. 
In the cylindrical case, the circumference of the cylinder 
provides the length scale that determines the compatibil- 
ity of candidate structures. 

In the next section, we evaluate the ground state prob- 
lem by considering a rather general set of alternative 
structures. In Section III, we compute the specific heat 
for several values of the ring number. Finally, in Sec- 
tion IV we estimate the quantum effects (by including 
zero-point energy) present in a rather extreme case, ^He 
atoms on a cylindrical surface, and compare the resulting 
energy with the known energy. 



GROUND STATE ENERGY 

We assume that the ground state structure is close- 
packed, as exemplified in Figure 2, and that it can be 
derived with the following algorithm. At any specific 
value, say zero, of the coordinate (z) parallel to the cylin- 
der's axis, there are v atoms, distributed uniformly in 
azimuthal angle 4>. We call v the ring number of the 
structure and consider structures with integral value of 
V. Figure 2 depicts the case v = A, which turns out to 
be an important example. A unit cell of this structure 
consists of four atoms at z=0 (at azimuthal angles 0=0, 
7r/2, TT and 37r/2) and four atoms at z=a (0 = 7r/4, 37r/4, 
57r/4 and 77r/4). The structure is characterized by a one- 
dimensional density p = 4/a, since there are 8 atoms in a 
unit cell of length 2a. For each such hypothetical struc- 
ture, we have performed two calculations. The first is a 
total energy calculation, aimed at determining the lowest 
energy structure. The energy in this case is taken as the 
sum of two-body L-J interactions between atoms at all 
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lattice sites. The second study yields the thermodynamic 
properties, described in the next section of this paper. 

Before embarking on these calculations, we assess the 
models. The only approximations in the ground state 
calculation are the use of Lennard- Jones potentials, the 
omission of kinetic energy and the assumption that the 
actual structure fits the close-packing description. The 
first two are conventional approximations in the lattice 
gas approach; we note only that many-body corrections 
may be important in this geometry (but we ignore them) 
The third assumption seems logical, since all sim- 
ple close-packed lattice structures are included. One can 
imagine other possibilities, such as one in which the unit 
cell consists of more than two rings, but these seem im- 
plausible to us. We note that the problem of packing on a 
spherical surface is quite different from that on a cylinder; 
there, frustration arises because of the difficulty of sat- 
isfying local packing requirements 0. Here, instead, 
we find many high density, strongly bound and nearly 
degenerate low energy structures, which do satisfy local 
bonding requirements. 

We have obtained the (ground state) energy results 
from the following procedure. For any assumed values of 
the ring number and cylinder radius, the energy E^{a, R) 
is evaluated by summing contributions from all inter- 
atomic interactions: 

E,{a,R)^Y.U{n,) (1) 

Here, the sum includes all pairs of atoms, separated by 
a 3D distance r^j. It is convenient to measure distances 
in units of the hard core parameter, to permit scaling 
between different solutions. Thus, R* and a* = a/ a are 
reduced distance variables; similarly E* = E/e is a re- 
duced energy. As an example (i^ — 4) revealing the en- 
ergy's dependence on these lengths, a contour plot of the 
function Ei{a* , R*) is shown in Figure 3. Note that this 
function possesses two local minima; these correspond to 
two quite distinct geometries. The minimum with the 
larger value of a* (smaller value of R*) corresponds to 
neighbors within the same ring that are separated by 
Ar ~ Tmin (the equilibrium distance of the pair poten- 
tial). The other minimum energy configuration (larger 
value of R*) involves nearest neighbors in adjacent rings, 
separated by Ar « rmi„ ||]. While, in either case, the 
low energy of the structure comes primarily from such 
optimization of nearest neighbor distances, longer range 
interactions do play a significant role in determining the 
total energy. This is evident from the fact that the (re- 
duced) cohesive energy per particle has a maximum value 
as high as 3.62 for ly ~ 4. This is 45 % greater than the 
value (5/2 for v ~ A) that would result if only nearest 
neighbor interactions were included. 

For each pair of values of v and R, one thus determines 
a unique value a*„j„ for which this energy function is a 



global minimum. This optimized value of a^j„ appears in 
the lower panel of Figure 4 and the corresponding energy 
appears in the upper panel. In the latter, one observes 
two alternative behaviors: either a single value of R yields 
a minimum in this function (for v—l or 2) or two values 
yield local minima (for v > 2). In the latter case, with 
the single exception of i/=3, the lower energy structure 
is the one with the smaller value of R, i.e., the case with 
nearest neighbors in the same ring. 

By such an analysis, we derive the ground state en- 
ergy, shown in Figure 5, representing the global lower 
bound of the ensemble of curves in Figure 4. This scal- 
lop shell-like function manifests the following plausible 
behavior. For very small R*, the lowest solution corre- 
sponds to the case v=l. Near R*=0.5, the lowest energy 
shifts to the v=2 structure; this is a plausible result be- 
cause then two atoms may occupy the same ring with- 
out significant hard core repulsion. For increasing R, the 
minimum energy and structure undergo a sequence of 
transitions between different values of v. Interestingly, 
the sequence is not monotonic: after a very narrow re- 
gion (0.60 < R* < 0.63) in which is optimal, there 
occurs a region (0.63 < R* < 0.73) in which 1^=2 is opti- 
mal, followed by an extended region (0.73 < R* < 0.87) 
in which i^=4 provides the lowest energy. Note that there 
are many energy minima close to the (reduced) energy 
E/{Ne) = —3.6. The global minimum energy structure 
occurs at R*=0.78, with reduced energy -3.62. At that 
point, there are two neighbors at reduced distance 1.10 
and four at distance 1.11, both of which are very close 
to the pair potential minimum value, r^i„ ~ 1.12. This 
most cohesive configuration corresponds to a cohesive en- 
ergy some 7 % higher than the two-dimensional ground 
state energy (3.38) of the L-J potential (corresponding 
to a hexagonal packing) 1^. This result implies that if 
atoms were to self-assemble on a surface of any shape, the 
cylindrical surface would be stable relative to the planar 
surface. While we have not performed the corresponding 
calculation for atoms on a spherical surface, we suspect 
that the energy in that case would be competitive with 
the present results Based on experience found in 

the case of a Coulomb interaction , we expect that 
frustration due to many energetically similar structures 
would be likely to occur in the spherical case. We note, 
for completeness, that the reduced cohesive energy in ID 
is slightly greater than one (1.03), while in 3D the value 
is quite large (8.1), a result of both the higher coordi- 
nation possible in 3D and the large contribution of long 
range forces. 

A structure of the type we are studying will sustain 
sound waves with various polarizations. The simplest 
such wave is a longitudinal compressional wave, with 
propagation vector parallel to the cylinder's axis. In the 
long wavelength limit, the corresponding sound speed for 
the case of mass M particles will satisfy: 
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Ms^ = , o (2) 

Here the derivative is evaluated at the ground state 
configuration for any R. Figure 6 shows corresponding en- 
ergy curves from which the derivatives in Equation 2 may 
be computed. For the curves shown, the reduced sound 
speed MsVe has the values 86, 37 and 91 at R*=0.78, 
1.05 and 1.25 (the three lowest energy minima in Figure 
6), respectively. In the case of a mass 16 particle (CH4), 
these values correspond to s w 3 km/s, comparable to 
the bulk speed of sound of CH4. The very high value is 
indicative of a very tightly bound and rigid structure. 

THERMODYNAMICS 

There exists a venerable tradition of applying the 
(Ising) lattice gas model to describe the condensation of 
gases. The critical exponents of the liquid- vapor tran- 
sition are believed to be exactly determined with this 
model. As is common in such applications, we simplify 
our calculations by assuming that only nearest neighbors 
interact and these interactions all have the same energy 
(—J). While J might simply be set equal to the well- 
depth of the pair potential, a more sophisticated model 
might increase the value of J to incorporate attractive, 
longer range interactions in some approximation. The 
previous section's results for the energy would suggest 
that an increase of J/e by a factor ~ 1.5 is needed to 
derive the ground state's energy. However, at a density 
less than complete filling of sites, the occupation fraction 
would reduce this hypothetical long range correction sig- 
nificantly. Of course, accurate calculations (not under- 
taken here) would incorporate longer range interactions 
explicitly in the Hamiltonian itself. 

Our method of study is the explicit evaluation of the 
partition function, within the canonical ensemble. For 
the case when N sites are occupied at temperature T, 
this function is 

Q^(T)=^e-^t"'>/^ (3) 

Here, Boltzmann's constant is taken to be one, the sum is 

over all configurations {rii} that yield a total of N parti- 
cles (out of Ns sites), and E{ni} is the corresponding en- 
ergy. Each configuration corresponds to a specific choice 
of occupied sites. Periodic boimdary conditions are em- 
ployed, so that the right end sites of a periodic region 
interact with "neighboring" left end sites of that region. 
A typical calculation involves a configuration determined 
by the occupancy of the 6v sites contained within three 
unit cells of the lattice. We explore the accuracy of this 
procedure by varying the size of the periodic cell. Be- 
cause there is no transition in this ID system, the finite 



size effects do not attenuate any divergence in the spe- 
cific heat, but they are observable. Checks on the results 
come from the known ground state energy and the high 
T energy, obtained from a random site occupancy: 

Here 7 is the coordination number and is the occu- 
pied fraction of sites. Another check comes from the en- 
tropy S{T), which is obtained by integrating the heat 
capacity, divided by T, from zero to infinity, where 
5 = ln{iV,!/[iV!(Ar,-iV)!]}. Note that S/N, = -ein9- 
(1 - 6*) ln(l - 0), i.e. In 2 at 50 % occupancy. 

Figures 7 to 9 show the energy and specific heat as 
functions of the reduced temperature T*=T/J for the 
cases v—2, 3 and 4, respectively. We make a number of 
remarks about the results. First, the curves are all qual- 
itatively similar; this is not surprising because none of 
these finite ly cases exhibits a phase transition. Hence, all 
of the interesting behavior is concentrated in the regime 
near and below T*=l. The key qualitative dependence is 
a concentration, with increasing i> of the thermal "activ- 
ity" into a progressively narrower region of T. This trend 
is plausible because the limit of very large u is the (tri- 
angular lattice) 2D limit, which exhibits a genuine phase 
transition, with critical temperature T* « 0.91. 

Size effects are present in our calculations and may 
also occur in nature, where nanotubes are finite or may 
have finite segments that are perfectly ordered. To ex- 
plore these effects, we focus on the case v = 3. Fig- 
ure 10 shows how the ground state energy (inset in left 
panel), at 0=1/2, depends on Ng, the number of sites 
in the unit cell, while the high T result remains fixed at 
E/{N J) = —3/2, the exact result at high T arising from 
the 50 % site occupancy. The dependence on Ng for 
61=1/2 conforms to the expression E{T = 0)/{-NJ) = 
3 — 12 /Ns, derived from the bulk and surface energies of 
a ground state "island" consisting of A^s/6 isolated rings. 
Hence, the largest system shown {Ng = 30) has a ground 
state energy differing from the infinite system result by 
a fraction 4/30 0.14. The specific heat bump moves 
to progressively higher T as Ng increases; it is seen to be 
converging to a well defined limit at that point, with a 
maximum near T*=0.89, as is expected from the limiting 
behavior described above. 

One particularly interesting feature is the presence of a 
double maximum in the specific heat at fractional occu- 
pations both near, but not equal to, 1/2 (Fig. 10, right 
panel) . Note that this behavior becomes increasingly ev- 
ident as the system size grows, indicating that it is not 
a finite size artifact. To explore this phenomenon, we 
compute two correlation functions, defined in the follow- 
ing way: The transverse correlation function is obtained 
from the average of the product of the occupation num- 
bers of the three sites in the same ring, and the longitu- 
dinal correlation function is the average of the product 
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of the occupation numbers of three sites in consecutive 
rings. 

These correlation functions are plotted in Figure 1 1 in 
the case of z/=3, for the same occupations as in Fig. 8. 
One observes drastically different behavior for — 1/2 
and 7^ 1/2. For 6 — 1/2, there arises transverse cor- 
relation below T*=0.75, persisting down to T=0, while 
the longitudinal correlations remains low and constant 
for the whole range of temperature. The ground state 
includes completely filled rings (perfect transverse order) 
but the longitudinal order is imperfect, e.g. because the 
islands have edges at 6* = 1/2. Note that the peak in 
the specific heat appears at the same temperature (T*« 
0.75) for which the transverse correlations start to de- 
velop, indicating a quasi transition to a more ordered 
state along the azimuthal direction. 

How do we explain the second, low T, bump in C for 6 
close to 1/2 ? We attribute it to development of the lon- 
gitudinal correlation at a lower T since the rings are not 
all complete in this case. The upper panels of Figure 11 
{9 — 8/18 and 10/18) show how the longitudinal corre- 
lation arises for T* below 0.25, i.e. a much lower T than 
where the transverse order begins, giving rise to the cor- 
responding low T bump in Figure 8. Rings start to form 
at T* « 0.75 and this ordering leads to a longitudinal 
quasi-condensation when these rings start to order along 
the z direction as T is lowered. When the occupation 
is considerably larger that 1/2 (for example as shown in 
the bottom right panel of Fig. 11, for = 15/18), the 
longitudinal order prevails because almost all the rings 
are occupied and the low T islands are bigger. 

Based on the behavior of these two correlation func- 
tions we can also explain the absence of two peaks in 
the V = 2 and v — A cases. For half filling, the longi- 
tudinal order is always poor because the system tends 
to aggregate in an island and the peak in the specific 
heat corresponds to the development of the azimuthal or- 
der characteristic of this structure. For occupancy much 
greater than 1/2, the system tries to condense along the 
z direction, which coincides with the formation of a much 
bigger island and the specific heat peak indicates the the 
appearence of this longitudinal order. However, because 
the energy difference between the transverse and longi- 
tudinally ordered structures is smaller than in the i/ = 3 
case, for intermediate cases when is close to 1/2, longi- 
tudinal and transverse ordering occur at nearly the same 
T, resulting in a single peak in the specific heat. 

QUANTUM EFFECTS 

The preceding discussion deals with the classical lattice 
problem, leaving open the question of quantum effects. 
If one were to compute the rms fluctuations of atomic 
displacements due to zero-point motion, the (essentially) 
ID phonon states of this lattice would lead to a diver- 



gent result. Hence, no crystalline state is possible for 
our lattice, even at T=0. However, as in the analogous 
2D problem at finite T, the divergence might not elimi- 
nate the possibility of a solid phase, exhibiting topolog- 
ical long range order (as distinct from crystalline order, 
i.e. periodicity) p], p^ . We do not speculate further 
about this problem here. 

However, we do wish to consider the implication of our 
results for the case of a rather extreme quantum problem, 
the ground state of ^He atoms confined to a cylindrical 
surface. In this case, there is no ambiguity about the 
ground state, since it is surely a liquid, as in 3D. Our 
concern here is a quantitative one: what is the binding 
energy of this ground state? We answer this question 
with an estimation method similar to one used by F. 
London more than half a century ago in describing ^He 
in 3D ||l^. The goal here is to see whether the lattice 
model contains any physics relevant to the *He case. 

One estimate of the total energy Etot of the system is 
the following: 

^^^1 =< V >o +i^est (5) 

Here < V >o is the ground state energy computed 
in Section II for a cylindrical lattice. The latter is a 
lower bound to the true potential energy of the quantum 
problem. Hence, the total iSg^] is a lower bound if the 
kinetic energy K^st were known to be less than, or equal 
to, the true kinetic energy of ^He. Unfortunately, we 
cannot establish such a relation without carrying out a 
full calculation. Instead, we approximate K^st, following 
London, leading to a total energy estimate, instead of a 
genuine bound. Our estimate of the kinetic energy per 
particle includes contributions from both azimuthal and 
longitudinal motions for atoms within the lattice: 

Figure 12 shows the results of this calculation, along 
with several others. One of these is a variational calcu- 
lation of the ground state energy Eyar of liquid ^He on 
a cylindrical surface, by Carraro |Q. While this latter 
is just an upper bound to the exact energy, it probably 
comes within 0.2 K of the exact energy. In comparing 
these curves, we note that both the energy estimates E^st 
and E^ar exhibit a common feature: there is a single to- 
tal energy minimum in the vicinity of R* « 0.7 to 0.8. 
We observe that iSg^] lies significantly below Carraro's 
results for all values of R*. The large discrepancy (w 15 
K) is attributable, in our opinion, primarily to the ne- 
glect of the large increase in potential energy above the 
classical lattice energy < V >q. To analyze this, we make 
a drastic approximation: that the motion of the system is 
totally harmonic, i.e. a phonon description. In that case, 
the potential energy increase above the value < 1/ >o is 
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equal to the kinetic energy per particle. This leads to a 
revised estimate of the energy: 

Ef}^ =< V >o +2ife.t (7) 

As seen in Figure 12, the variational results lie midway 
between Eg^2 ^^'^ ^est ^^'^ minima of these functions 
lie close to that of the variational calculations. This com- 
parison suggests that our model is getting the physics 
approximately correct. However, one must bear in mind 
that a more complete analysis would take into account 
the relaxation of the system to incorporate the kinetic 
energy. This would be the analog of allowing the liquid 
density to vary in searching for the ground state, as is 
conventionally done in variational calculations (and was 
done by London with his analogous model). The result 
of such an approach would be a reduction in the density 
(as occurs for 3D ^He and H2) and in the magnitudes 
of all energies in the system. Such an extended analysis 
seems unwarranted in view of the naivete of the present 
description and the availability of alternative, more ac- 
curate, computational methods. 

We are grateful to Carlo Carraro for helpful discussion. 
This research has been supported by the Petroleum Re- 
search Fund of the American Chemical Society and by 
Fundacion Antorchas. 
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FIG. 1: Density of H2 molecules at T=10 K as a function of radial distance inside a nanotube of radius 7 A for several values 
of the chemical potential /u (taken from Ref. [2]). 
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FIG. 2: Schematic depiction of the cylindrical lattice structure having ring number v=4. The lattice constant is 2a. 
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FIG. 4: Upper panel shows the energy per particle (in units of the pair potential's well depth) as a function of reduced radius, 
for various assumed ring numbers, v—l, 2, 3, 4, 5, 6, and 8, from left to right. Each value of E is derived by choosing the 
optimized lattice constant {a,„i„), shown in the lower panel. The odd v curves are dashed and v—l, 2, 6 and 8 are labelled. 




FIG. 5: Energy as a function of radius, obtained by selecting the minimum energy from the solutions in Figure 4. The global 
minimum energy for this problem occurs at R*=0.78 = 4), as indicated by the happy face. 




FIG. 6: Energy as a function of a* for various values of R* in the case v = 4, see Fig. 3. The curve corresponding to R*=0.78 
yields the global energy minimum of the problem. 




FIG. 7: Energy (lower panel) and specific heat (upper panel) for the case of ring number v = 2, Ng = 16 (8 rings). The 
full curve corresponds to half filling occupancy (0=8/16), the dotted curve to 6=7/16, the dashed one to 0=9/16, and the 
dot-dashed to 61=13/16. 




FIG. 8: Energy (lower panel) and specific heat (upper panel)for the case of ring number v = 3, Ns = 18 (6 rings). The full curve 
corresponds to half filling occupancy (^=9/18), the dotted curve to 0=8/18, the dashed one to 0=10/18, and the dot-dashed 
to 6'=15/18. 
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FIG. 9: Energy (lower panel) and specific heat for the case of ring number v = 4, Ns = 16 (4 rings). The full curve corresponds 
to half filling occupancy (^=8/16), the dotted curve to 0=7/16, the dashed one to ^=9/16, and the dot-dashed to 0=13/16. 
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FIG. 10: Specific heat for v — 3, showing the size effect as the number of sites Ns increases. The full curve depicts the case 
for A'^s = 18 (6 rings), the dashed curve corresponds to Ns = 24 (8 rings) and the dotted curve to Ns =30 (10 rings). The 
left panel show the case for half occupancy (the inset depicts the energy per particle) whereas the right panel corresponds to 
occupancy below 1/2: 6 = 0.44 (full curve), 0.42 (dashed curve), 0.43 (dotted curve). 




FIG. 11: Transverse (full curves) and longitudinal (dashed curves) correlation functions (defined in the text) for the same 
occupancies as in Fig. 8 {v = 3). 




R/a 

FIG. 12: The ground state energy of ''He atoms confined to a cylindrical surface. The full curve shows the result of the 
variational calculation i?«ar by Carraro (Ref. 1). Also shown are two alternative estimates E^]^l and E^^l of the total energy 
obtained from the classical potential energy )V(o and estimated kinetic energy, K^st, as described in the text. 



